Molecular Dynamics Study of a Thermal Expansion Coefficient: 
Ti Bulk with an Elastic Minimum Image Method 
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Linear thermal expansion coefficient (TEC) of Ti bulk is investigated by means of molecular 
dynamics simulation. Elastic Minimum Image Convention of periodic boundary conditions is intro- 
duced to allow the bulk to adjust its size according to the new fixed temperature. The TEC and 
the specific heat of Ti are compared to the available theoretical and experimental data. 
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INTRODUCTION 



ELASTIC MINIMUM IMAGE CONVENTION 
(EMIC) 
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Titanium is not present in pure form in the nature, 
and its composites have different thermal expansion co- 
efficients (TEC) Q. Determining the TEC values may 
help to guess the composition. For expansion simula- 
tions, various Molecular Dynamics (MD) methods allow- 
ing for a volume change exist 

Hi- 

Expansion is defined 
through a change of the bond length which may be de- 
termined with several methods, like the average length 
of the bulk diagonals For a more statistically de- 
pendable result, the mean distances between the atoms 
in the bulk are calculated in this work. Expansion of 
the sample requires elasticity of the periodic boundary 
conditions (PBC), but the PBC and its minimum image 
convention use fixed box lengths in the simulation [2|. 
Therefore, elastic minimum image convention (EMIC) of 
the PBC is introduced. It resembles Berendsen's baro- 
stat method. However, Berendsen's method allows only 
isotropic changes in the volume of the simulation box [2J , 
while the EMIC method allows change of the shape as 
well as its size. For the present purpose correctness of the 
elastic constants and phonon frequencies are important. 
The recent many-body potentials such as Finnis-Sinclair, 
Gupta, and Glue potentials and embedded atom method 
cannot reproduce all elastic constants correctly [B| , while 
Lennard- Jones potential (LJ) Q give good results Q for 
Ti. Hence the LJ potential is employed into our tem- 
perature scaled isoenergetic and isobaric MD codes 
Relaxation runs of the bulk Ti have been performed by 
using this code. The Verlet integration algorithm Q is 
used, since it allows for time- reversible actions Q. 
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The EMIC method is used to adjust the simulation 
box lengths to the temperature changes during a run [J] . 
The entire simulation time ttot is divided into nb number 
of time-blocks Tb, each consisting of m time-steps ts, 
with St time step size (tss), which can be expressed as 
Uot — nb x Tb and Tb = m x St. 

The mean bond length M BL is calculated from 
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where N, C and are the total number of atoms, com- 
bination of N by 2 and pair distance in the simulation 
cell, respectively. The changes in MBL are reflected di- 
rectly to the box length. Although, temperature scaling 
is applied after each Tb, EMIC is only applied after some 
number oiTb, let's call nscb. Therefore, EMIC is applied 
at every ecb time interval, i.e., ecb — nscb x Tb, where 
nscb = [1-100]. The system temperature fluctuates freely 
during each Tb, while the MBL differences accumulate 
and the system tries to reach an equilibrium during each 
ecb [4| . In this way, the system also fulfils ergodicity. The 
difference in MBL during the k th and (k — l) th ecb is 
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where MBLi is the average during the i Tb, and the 
terms on the right hand side are the MBL averages at 
the k th and (k - l) th ecb. An effective MBL average at 
the k th ecb is obtained by multiplication of AMBL kt k-i 
with an adjustable parameter p as 
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where < p < 1. Note that Eq. [3] becomes Eq. [2] 
if p = 1. Finally, the box length of the EMIC in the 
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x-direction from k to (fc + 1) ecb is 
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Similar equations hold for the other directions. Since the 
radial MBL averages are reflected to each of the direc- 
tions as in Eq. [31 the initial shape of the bulk is forced 
remain to unchanged. However, the real lattice may not 
expand with the same rate along each of the Cartesian 
directions 13, II | ■ Therefore Eq. [3] is reformulated for 
each coordinate separately as 

(MB^r/Jik) =P&MBLx h , k - X + (MBm)l{l {k _ 1) (5) 

Similarly, Eq. [5] is also reformulated for each direction 
separately as 
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Values of m in a Tb should be chosen as large as possible 
for the numerical convergence of the physical quantities, 
in order to reduce non-physical effects [2j . Hence, several 
test runs are required to determine the parameters m, 
nscb, and p. 



SIMULATION 

The conventional a-Ti has hep structure. The charac- 
teristics of a hep cell are determined by proportionality 
of simple cell sizes c/a. This ratio and the binding energy 
parameter £o are taken as 1.62073 and 0.20154 eV when 
producing the L J potential for Ti Q , while the recent ex- 
perimental c/a value is 1.601 12]. The unit cell sizes are 
a = 2.94055 A and c = 4.76583 A for LJ potential com- 
pared to a = 2.953 A and c = 4.729 A from experiment. 
The bulk sample in the simulation is prepared using the 
former set of constants. A safe cut-off radius for the inter- 
actions is chosen as 2.8 a (a = 2 1 / 6 r , o) 0, where the bond 
length r is 2.96292 A [7]. Since the box length should 
be equal or greater than twice the cut-off radius , the 
smallest possible sample having 389 Ti atoms with a lat- 
tice size of 6a x 66 x 4c, i.e., 17.718 x 17.718 x 18.916 A 3 
is constructed. 

Simulations have been carried out for the temperature 
range 100-400 K, in steps of 100 K. Quantum effects on 
the specific heat become important at low temperatures 
[l3j |; as these effects are not included in the binary po- 
14j . the lower temperature was chosen as 100 



tentials 

K. The higher temperature limit is set near the Debye's 
temperature, since the binary potentials require modifi- 
cations for the higher temperatures Q ■ The final veloci- 
ties and coordinates of the previous run are used as the 
initial conditions for the next run, i.e., those of the 200 



K are used as the initial values for the 300 K run, and so 
on. At each temperature two different simulations were 
performed one after the other. Firstly, the bulk is run 
under constant pressure by using the EMIC boundary 
conditions to adjust the bulk sides. After that, the code 
is run by using PBC under this constant volume. 

The time averaged temperature over each Tb is calcu- 
lated from 
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where (E k i„) is the time averaged kinetic energy of the 
atoms, N is the number of atoms in the sample, and 
ks is the Boltzmann constant. Hence, the simulation 
temperature is the average over all T(K) values. 

A number of runs have been made to find the right 
parameters. The tss and Tb are varied between 0.2-2 fs 
(lfs = 10~ 15 sec) and 10-100 tss, respectively. The total 
number of time steps lie within the range of 163,000- 
2,002,000 and the resulting total run times are varied in 
between 18,700 fs and 402,000 fs. The coefficient p of 
the effective MBL differences is used as 0.5-1 in order to 
let the system change smoothly. Otherwise, if p is taken 
equal to 1 , a sudden change in the periodic box size may 
lead to excess energy at some particles. 

Firstly, average bond differences AMBL kk -i are dis- 
tributed proportionally, according to c/a ratio, to each 
of the Cartesian directions of the EMIC box by Eqs. [3] 
and [4j Therefore the shape of the bulk is enforced to 
remain globally unchanged. However, this brutal shar- 
ing brought the problem of exceeding force accumulation, 
and the problem of stability at 300 K and 400 K. Using 
small time step sizes has solved the problem of 300 K, 
but not at the higher temperatures. Then, AMBL kk _i 
is shared to each Cartesian coordinate separately, and 
the bulk is reconstructed in each direction independently 
(Eqs. E]and[|5]). This causes accidental image conventions 
of the coordinates after some number of time intervals. 
In order to avoid these and to reach numerical conver- 
gence of the physical quantities, some number of runs 
have been done to determine the run parameters. 

Together with the temperature of the bulk, coordinates 
of the atoms along each of the Cartesian directions, ki- 
netic, potential and total energies of the atoms are calcu- 
lated. Fluctuations on all these quantities, specific heat 
and thermal expansion coefficients are obtained. 

The linear thermal expansion coefficient (TEC) is de- 
fined as 
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where Lq, AL, and AT are the initial length, change of 
length and change of temperature, respectively. In prac- 
tice, L is calculated with various methods throughout the 
entire run. In this work, MBL and MBLq are used in 
place of L and Lq, respectively. 
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Heat capacity can be observed by doing a series of 
simulations at different temperatures with correspond- 
ing constant volumes. It can be calculated in each run 
separately by using temperature or energy fluctuations. 
Various authors have offered similar forms for the calcu- 
lation of heat capacity 0, [§] . It is defined through kinetic 
energy fluctuations for the microcanonical ensembles by 
Lebowitz et al. [la as 
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where a 2 (KE) is mean square deviation of the kinetic 
energy throughout the entire run. 



RESULTS AND DISCUSSIONS 
Shape of the Bulk 

Small changes on the coordinates of the bulk atoms 
are observed. Typical bond length changes are about 
0.01-0.04 A (0.34-1.35% of T = K bond length). The 
centre of mass is unchanged during the simulations. Be- 
tween 100-300 K the shape of the bulk remains more or 
less unchanged. However, at 400 K the resulting bulk has 
a shape anomaly like in the experiments where the ex- 
pansion of the bulk varies in different directions fiol . 11 1 ; 
the difference of the c/a value after the 400 K simulation 
is roughly 0.4% in comparison to its T = K value. 
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FIG. 1: Caloric curve of the Ti-389 bulk obtained by the 
simulation. 



0. Moreover, the temperature values are observed to be 
certain within [0.004, 0.047] K error. The heat capacities 
are observed by using the method explained in the pre- 
ceding section. These are calculated at each temperature 
using Eq. [5] The results are compared to the available 
experimental data H, 17| in Fig. [5J There exists a good 
agreement. The Maximum difference occurs at 100 K as 
of 17%. The mean square deviation of the kinetic energy 
of the atoms in Eq. (11) is taken as Is, if it were taken as 
2s which covers the 95% of the region- they would have 
certainly covered the whole experimental values. 
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FIG. 2: Comparison of this simulation and experimental re- 
sults for the specific heat capacities of Ti. • present simula- 
tion; A, experiment, Ref. [lg|; V, experiment, Ref. [l7| . 
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FIG. 3: The calculated rate of change of the mixed bond 
length values against temperature. The solid curve represents 
the linear best fit. 



Specific Heat Capacity 

The temperature development at different total ener- 
gies, resulted from the simulation, is shown in Fig. [TJ A 
linear behaviour and a small increase are seen by increas- 
ing the temperature, although the total energies at each 
temperature are constant with small fluctuations. The 
maximum fluctuations are roughly 10 -5 eV as expected 



Thermal Expansion Coefficient 

The linear thermal expansion coefficient of the Ti is 
calculated from Eq. [8] using the MBL. The MBL varies 
with the temperature of the bulk. The temperature de- 
pendence of the relative MBL change, AMBL/MBL , 
is shown in Fig. [3J where MBL$ is the mixed bond length 
at T = K. It is important to note that the T = K 
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value is only a reference point and not a simulation re- 
sult. As seen, the curve exhibits a nearly linear increase 
in this temperature range. 



TABLE I: Experimental and our simulation results (in units 
of 10~ 6 /K) for the linear thermal expansion coefficient of Ti 
as function of temperature T. 



T (K) 


MD 


Experiments 


too 


4.48 




200 


7.73 








8.35 (a) 






8.5 (b) 


300 


8.50 








8.6 (c) 






11.9 (d) 


400 


10.99 


9.32 (b) 


300-400 


9.74 (f) 


8.7 (c) 


(a) T= 298.f5 K, Ref. [18] 


(b) Ref. [11] 




(c) T= 298.15 K, Ref. [12] 


(d) T= 23 °C (296.15 K), Ref. [19] 


(e) T= 20-100 °C (293.15-393.15 K), Ref. [18] 


(f) Mean value of simulations between T — 300-400 K. 



The calculated TEC of the Ti values are tabulated in 
Table U together with the available experimental results. 
We could not find any experimental TEC value below 
300 K to compare with our simulation results. The ex- 
perimental results at 300 K [ll], [lj, El, Eil show notice- 
able irregularity One of them perfectly matches with 
the simulation result. At 400 K, the simulation result of 
10.99 x 10 _6 /K differs by 18% from the experimental re- 



sult of Ref. [11|]. The difference is smaller in the 300-400 
K range, 12%, where the simulation and experimental 
results are 9.74 x 10~ 6 /K and 8.7 x 10~ 6 /K, respectively. 



SUMMARY 

Isobaric and isoenergetic MD simulations have been 
performed with Verlet algorithm. The LJ potential pa- 
rameters of Yamamoto-Kagawa-Doyama is used for the 
Ti-Ti interactions. A Ti-389 sample is produced by us- 
ing the same lattice constants of this potential. The elas- 
tic minimum image convention of the periodic boundary 
conditions is set to enable the elasticity of the bulk. Con- 
stant pressure, and isoenergetic simulations have been 
done, one after the other. Specific heat and linear TEC 
of Ti are calculated between 100 K and 400 K. The lower 
temperature is determined to reduce the quantum effects 
[HI]. And the upper limit is set to the Debye's tempera- 
ture since the binary potentials require modifications for 
the higher temperatures The calculated TEC value 
at 300 K matches perfectly with one of the experiment, 



while the highest difference occurred as 18% at 400 K. 
Furthermore, observed path of the specific heat in the 
simulation has a similar path with the experiments as 
shown in Fig. [21 Numerical convergence problems at 
both temperature limits are observed in the simulation. 
Solving the high temperature problems will be our future 
work to clarify the behaviour in this regime. 
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